Exploration des données de vent RDPS

import xarray as xr
import xesmf as xe
import hvplot.xarray
import pandas as pd
import os
import glob
import geopandas as gpd
from distributed import Client
import matplotlib.pyplot as plt

os.chdir('/Users/caramelo/Documents/GitHub/cmc/')
grib_list_total=sorted(glob.glob('test_reg*.grib2'))
#grib_list=grib_list_total[0:7]
#grib_list

grib_list_total

['test_reg202103260.grib2',
 'test_reg202103261.grib2',
 'test_reg2021032610.grib2',
 'test_reg2021032611.grib2',
 'test_reg2021032612.grib2',
 'test_reg2021032613.grib2',
 'test_reg2021032614.grib2',
 'test_reg2021032615.grib2',
 'test_reg2021032616.grib2',
 'test_reg2021032617.grib2',
 'test_reg2021032618.grib2',
 'test_reg2021032619.grib2',
 'test_reg202103262.grib2',
 'test_reg2021032620.grib2',
 'test_reg2021032621.grib2',
 'test_reg2021032622.grib2',
 'test_reg2021032623.grib2',
 'test_reg2021032624.grib2',
 'test_reg2021032625.grib2',
 'test_reg2021032626.grib2',
 'test_reg2021032627.grib2',
 'test_reg2021032628.grib2',
 'test_reg2021032629.grib2',
 'test_reg202103263.grib2',
 'test_reg2021032630.grib2',
 'test_reg2021032631.grib2',
 'test_reg2021032632.grib2',
 'test_reg2021032633.grib2',
 'test_reg2021032634.grib2',
 'test_reg2021032635.grib2',
 'test_reg2021032636.grib2',
 'test_reg2021032637.grib2',
 'test_reg2021032638.grib2',
 'test_reg2021032639.grib2',
 'test_reg202103264.grib2',
 'test_reg2021032640.grib2',
 'test_reg2021032641.grib2',
 'test_reg2021032642.grib2',
 'test_reg2021032643.grib2',
 'test_reg2021032644.grib2',
 'test_reg2021032645.grib2',
 'test_reg2021032646.grib2',
 'test_reg2021032647.grib2',
 'test_reg2021032648.grib2',
 'test_reg2021032649.grib2',
 'test_reg202103265.grib2',
 'test_reg2021032650.grib2',
 'test_reg2021032651.grib2',
 'test_reg2021032652.grib2',
 'test_reg2021032653.grib2',
 'test_reg2021032654.grib2',
 'test_reg2021032655.grib2',
 'test_reg2021032656.grib2',
 'test_reg2021032657.grib2',
 'test_reg2021032658.grib2',
 'test_reg2021032659.grib2',
 'test_reg202103266.grib2',
 'test_reg2021032660.grib2',
 'test_reg2021032661.grib2',
 'test_reg2021032662.grib2',
 'test_reg2021032663.grib2',
 'test_reg2021032664.grib2',
 'test_reg2021032665.grib2',
 'test_reg2021032666.grib2',
 'test_reg2021032667.grib2',
 'test_reg2021032668.grib2',
 'test_reg2021032669.grib2',
 'test_reg202103267.grib2',
 'test_reg2021032670.grib2',
 'test_reg2021032671.grib2',
 'test_reg2021032672.grib2',
 'test_reg2021032673.grib2',
 'test_reg2021032674.grib2',
 'test_reg2021032675.grib2',
 'test_reg2021032676.grib2',
 'test_reg2021032677.grib2',
 'test_reg2021032678.grib2',
 'test_reg2021032679.grib2',
 'test_reg202103268.grib2',
 'test_reg2021032680.grib2',
 'test_reg2021032681.grib2',
 'test_reg2021032682.grib2',
 'test_reg2021032683.grib2',
 'test_reg2021032684.grib2',
 'test_reg202103269.grib2']

ds=xr.open_mfdataset(grib_list_total,concat_dim='valid_time',engine='cfgrib',combine='nested')

ds

<xarray.Dataset>
Dimensions:            (valid_time: 85, x: 935, y: 824)
Coordinates:
    time               datetime64[ns] 2021-03-26
    step               (valid_time) timedelta64[ns] 00:00:00 ... 0 days 09:00:00
    heightAboveGround  int64 10
    latitude           (y, x) float64 dask.array<chunksize=(824, 935), meta=np.ndarray>
    longitude          (y, x) float64 dask.array<chunksize=(824, 935), meta=np.ndarray>
  * valid_time         (valid_time) datetime64[ns] 2021-03-26 ... 2021-03-26T...
Dimensions without coordinates: x, y
Data variables:
    gust               (valid_time, y, x) float32 dask.array<chunksize=(1, 824, 935), meta=np.ndarray>
Attributes:
    GRIB_edition:            2
    GRIB_centre:             cwao
    GRIB_centreDescription:  Canadian Meteorological Service - Montreal 
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             Canadian Meteorological Service - Montreal 
    history:                 2021-03-26T13:58:15 GRIB to CDM+CF via cfgrib-0....

ds.gust.attrs

{'GRIB_paramId': 260065,
 'GRIB_shortName': 'gust',
 'GRIB_units': 'm s**-1',
 'GRIB_name': 'Wind speed (gust)',
 'GRIB_cfVarName': 'gust',
 'GRIB_dataType': 'af',
 'GRIB_missingValue': 9999,
 'GRIB_numberOfPoints': 770440,
 'GRIB_typeOfLevel': 'heightAboveGround',
 'GRIB_NV': 0,
 'GRIB_stepUnits': 1,
 'GRIB_stepType': 'instant',
 'GRIB_gridType': 'polar_stereographic',
 'GRIB_gridDefinitionDescription': 'Polar stereographic ',
 'long_name': 'Wind speed (gust)',
 'units': 'm s**-1'}

ds.gust.isel(valid_time=1).plot()

<matplotlib.collections.QuadMesh at 0x167e8c9a0>

ds = ds.rename({'latitude': 'lat', 'longitude': 'lon'})
#ds
ds_out = xe.util.grid_global(0.1, 0.1)
#ds_out  # contains lat/lon values of cell centers and boundaries.
regridder = xe.Regridder(ds, ds_out, 'bilinear')
regridder
dr_out = regridder(ds)
dr_out.where(dr_out>0).gust.plot()
dr_out['lat'] = dr_out.lat[:,0].drop('lon')
dr_out['lon'] = dr_out.lon[0,:].drop('lat')
dr_out = dr_out.sortby(['x','y'])
dr_out = dr_out.rename({'x':'longitude',
                        'y':'latitude',
                        'lon':'longitude',
                        'lat':'latitude'})
dr_out = dr_out.where(dr_out>0) # remplacer 0 par nan

/Users/caramelo/anaconda3/envs/test_xesmf/lib/python3.9/site-packages/xesmf/util.py:94: UserWarning: 360 cannot be divided by d_lon = 0.1, might not cover the globe uniformally
  warnings.warn(
/Users/caramelo/anaconda3/envs/test_xesmf/lib/python3.9/site-packages/xesmf/util.py:100: UserWarning: 180 cannot be divided by d_lat = 0.1, might not cover the globe uniformally
  warnings.warn(
/Users/caramelo/anaconda3/envs/test_xesmf/lib/python3.9/site-packages/xesmf/frontend.py:487: FutureWarning: ``output_sizes`` should be given in the ``dask_gufunc_kwargs`` parameter. It will be removed as direct parameter in a future version.
  ds_out = xr.apply_ufunc(

dr_out

<xarray.Dataset>
Dimensions:            (latitude: 1800, longitude: 3600, valid_time: 85)
Coordinates:
    time               datetime64[ns] 2021-03-26
    step               (valid_time) timedelta64[ns] 00:00:00 ... 0 days 09:00:00
    heightAboveGround  int64 10
  * valid_time         (valid_time) datetime64[ns] 2021-03-26 ... 2021-03-26T...
  * longitude          (longitude) float64 -179.9 -179.9 -179.8 ... 179.8 179.9
  * latitude           (latitude) float64 -89.95 -89.85 -89.75 ... 89.85 89.95
Data variables:
    gust               (valid_time, latitude, longitude) float64 dask.array<chunksize=(1, 1800, 3600), meta=np.ndarray>
Attributes:
    regrid_method:  bilinear

os.chdir('/Users/caramelo/Documents/GitHub/cmc/')


coords_stations=pd.read_csv('Stations matrice scribe.csv')
coords_stations

	OACI	NOMSTN	LAT	LON
0	CYUL	MONTREAL (DORVAL)	45.4667	-73.7500
1	CYOW	OTTAWA	45.3167	-75.6667
2	CYVO	VAL D'OR	48.0500	-77.7833
3	CYSC	SHERBROOKE	45.4333	-71.6833
4	CYQB	QUEBEC	46.8000	-71.4000
5	CYRJ	ROBERVAL	48.5167	-72.2667
6	CYGL	LA GRANDE RIVIERE	53.6333	-77.7000
7	CYBC	BAIE COMEAU	49.1333	-68.2000
8	CYZV	SEPT-ILES	50.2167	-66.2500
9	CYRQ	TROIS-RIVIÈRES (A)	46.3539	-72.6792
10	CYMX	MIRABEL	45.6833	-74.0333
11	CWPK	PARENT	47.9167	-74.6167
12	CYYY	MONT JOLI	48.6167	-68.2167
13	CYWK	WABUSH LAKE	52.9333	-66.8667
14	CYAH	LA GRANDE IV	53.7500	-73.6667
15	CYYB	NORTH BAY	46.3667	-79.4167
16	CYGP	GASPE	48.7833	-64.4833

#Montréal
gust1d = dr_out.gust.sel(latitude=coords_stations.iloc[0,2], longitude=coords_stations.iloc[0,3],method='nearest')
df = gust1d.reset_coords(drop=True).to_dataframe()
df.plot()
plt.title(coords_stations.iloc[0,1])

Text(0.5, 1.0, 'MONTREAL (DORVAL)')

#bon maintenant admettons pour la station de sherbrooke
gust1d = dr_out.gust.sel(latitude=coords_stations.iloc[3,2], longitude=coords_stations.iloc[3,3],method='nearest')
df = gust1d.reset_coords(drop=True).to_dataframe()
df.plot()
plt.title(coords_stations.iloc[3,1])

Text(0.5, 1.0, 'SHERBROOKE')

Test avec Holoviews

import holoviews as hv
from holoviews import opts
hv.extension('bokeh')

import hvplot.pandas  # noqa
df.hvplot(title='Sherbrooke')

import hvplot.pandas  # noqa
df.hvplot(title='Montréal')